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ABSTRACT 


A  finite  element  formulation  with  linear  triangular  elements  was  used  to 
solve  the  steady-state,  two-dimensional  conduction  heat  transfer  equation 
in  the  condenser  wall  section  of  an  internally  finned  rotating  heat  pipe. 

A  FORTRAN  program  using  this  method  was  coupled  with  the  ADS  program 
for  automated  design  of  the  internal  heat  pipe  fin  geometry  to  optimize 
heat  transfer.  An  increase  in  surface  area,  which  increases  heat  transfer, 
also  increases  the  condensate  level,  which  decreases  heat  transfer. 

The  additional  condensate  level  does  not  offset  the  advantage  gained  by 
the  increased  surface  area.  The  investigation  provided  combinations  of 
fin  half  angle,  number  of  fins,  and  fin  height  for  an  optimum  design. 

Water  is  used  as  the  working  fluid  and  the  heat  pipe  is  constructed  from 
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I.  niTRODOCTION 


A.  THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a  closed  container  designed  to  transfer  a 
large  amount  of  heat  in  rotating  machinery.  Since  the  heat  pipe  operates 
on  a  closed  two-phase  cycle,  the  heat  transfer  capacity  is  greater  than  for 
solid  conductors.  Also,  the  thermal  response  time  is  less  than  with  solid 
conductors.  The  three  major  elemental  parts  of  the  rotating  heat  pipe  are: 
a  cylindrical  evaporator,  a  truncated  cone  condenser,  and  a  fixed  amount 
of  working  fluid  as  shown  in  figure  1. 

An  annulus  is  formed  by  the  working  fluid  in  the  evaporator.  This 
occurs  at  rotationary  speeds  above  the  critical  speed.  The  additicn  of  heat 
to  the*  evaporator  vaporizes  the  working  fluid.  A  pressure  differential 
between  the  evaporator  and  the  condenser  causes  the  vapor  to  flow 
towards  the  condenser.  The  vapor  is  transported  to  the  condenser  with  its 
latent  heat  of  vaporization.  Condensation  of  the  vapor  on  the  inner  wall  is 
caused  by  external  cooling.  This  condensation  releases  the  latent  heat  of 
evaporation.  This  condensate  is  forced  to  flow  back  to  the  evaporator  by 
a  component,  acting  along  the  condenser  W2dl,  of  the  centrifugal  force 
which  is  caused  by  the  rotation  of  the  heat  pipe.  As  the  condensate 
collects  in  the  evaporator  the  cycle  is  repeated. 

Since  the  evaporator  and  condenser  portions  of  a  heat  pipe  function 
independently,  needing  only  common  liquid  and  vapor  streams,  the  area 
over  which  heat  is  introduced  can  differ  in  size  and  shape  from  the  area 
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over  which  it  is  rejected,  provided  that  the  rate  at  which  the  liquid  is 
vaporized  does  not  exceed  the  rate  at  which  it  can  be  condensed. 
Therefore,  high  heat  fluxes  generated  over  relatively  small  areas  can  be 
dissipated  over  larger  areas  with  reduced  heat  fluxes;  allowing  a 
cylindriced  evaporator  and  a  truncated  cone  condenser. 

Capillaury  action  acts  to  drive  the  condensate  back  to  the  evaporator 
in  a  conventional  heat  pipe.  No  Kmitation  due  to  capillary  action  is 
encountered  in  a  rotating  heat  pipe  nor  are  externad.  pumps  or  gravity 
depended  on  for  the  flow  of  the  working  fluid.  Therefore,  the  rotating  heat 
pipe  can  be  used  in  any  orientation  [Ref.  1]. 

B.  OPERATING  LIMITS  OP  A  ROTATING  HEAT  PIPE 

The  first  theoretical  investigation  of  the  rotating  heat  pipe  conducted 
at  the  Naval  Postgraduate  School  was  performed  by  Ballback  [Ref.  2]  in 
1969.  Various  fluid  dynamic  mechanisms  limit  the  performance  of  a  rotating 
heat  pipe.  Ballback  [Ref.  2]  studied  these  meclianisms  and  an  estimation  of 
the  sonic  limit,  boiling  limit,  entrainment  limit,  and  condensing  limit  of 
performance  was  made. 

1.  The  Sonic  Limit 

The  maximum  flow  of  the  vapor  is  set  by  the  choked  flow 
condition  in  the  rotating  heat  pipe.  This  limiting  vapor  flow  rate  occurs 
when  the  heat  flux  is  increased  and  limits  the  amount  of  energy  the  vapor 
can  t'"»nsport.  The  rotating  heat  pipe  effectiveness  is  reduced  due  to  this 
limitation.  The  limiting  heat  transfer  rate  due  to  this  condition  is: 
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(1) 


<?»  - 

and  the  vapor  velocity  is  considered  to  be  sonic, 

U,-c-  (gJcRD^ 

where 

Uy  =  velocity  of  the  vapor  in  ft/sec,  and 
A  =  cross  sectional  area  for  the  vapor  flow  in  ft 
c  =  sonic  velocity  in  ft/sec 
gg  =  gravitational  constant 
k  =  ratio  of  specific  heats 
R  ~  gas  constant  in  ft-lbf/lbm  R 
T  *  absolute  temperature  in  degrees  Rankine 
py  =  density  of  vapor  in  Ibm/ft^ 

2.  The  Boiling  Limit 

The  transition  from  nucleate  to  film  boiling  was  hypothesized  by 
Kutateladze  [Ref.  3]  to  be  a  completely  hydrodynamic  process.  He 
determined  the  following  theoretical  formula  for  predicting  the  burnout 
flux: 


<?,  -  j-f  ;))* 

where 

K  =  constant  value 

y 

A];  =  heat  transfer  area  in  the  boiler  in  ft 
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hfg  =  latent  heat  vaporization  in  Btu/lbm 
o  =  surface  tension  in  Ibf/ft 
q  -  acceleration  of  gravity  in  ft/hr"^ 
pf  =  density  of  fluid  in  Ibm/ft^ 

Pv  =  density  of  vapor  in  Ibm/ft  . 

A  constant  value  for  K  in  the  range  of  0.13  to  0.19  was  suggested  by  the 
experimental  data  obtained  by  Kutateladze. 

3.  The  Entrainment  Lindt 

The  flooding  constraint  in  a  wickless  heat  pipe  was  examined  by 
Sakhuja  [Ref.  4].  The  correlation  he  developed  is: 


V,  ■  - - ^ - J -  (4) 

where 

s  heat  transfer  rate  in  Btu/hr 
Aj  =  flow  rate  in  ft^ 

C  =  dimensionless  constant,  0.725  for  tube  with  sharp  edged  flange 
hfg  =  latent  heat  of  vaporization  in  Btu/lbm 
g  =  acceleration  due  to  gravity  in  ft/hr 
D  =  inside  diameter  of  heat  pipe  in  ft 
pf  =  density  of  the  fluid  in  Ibm/ft 
py  =  density  of  the  vapor  in  Ibm/ft^. 
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4.  The  Condensing  Limit 

The  condenser  section  of  a  rotating  heat  pipe  was  modeled  as  a 
truncated  cone  by  Ballback  [Ref.  2].  Using  this  model,  the  condensation 
limitation  for  a  rotating  heat  pipe  was  determined  by  Ballback  [Ref.  2].  He 
developed  the  following  condensation  limit: 

where 

Qf  s  total  heat  transfer  rate  in  Btu/hr 

kf  =  thermal  conductivity  of  the  condensate  film  in  Btu/hr-ft-F 
pf  ~  density  of  fluid  in  Ibm/ft^ 

0)  s  angular  velocity  in  1/hr 

hfg  ~  latent  heat  of  vaporization  in  Btu/lbm 

Ts  =  saturation  temperature  in  F 

Tj,  =  inside  wall  temperature  in  F 

Uf  =  viscosity  of  fluid  in  Ibm/ft-hr 

^  =  half  cone  angle  in  degrees 

Rg  =  minimum  wall  radius  in  ft 

L  =  length  along  the  wall  of  the  condenser  in  feet. 

The  geometry  and  speed  of  the  rotating  heat  pipe,  and  the  physical 
properties  of  the  working  fluid  comprise  the  condensing  limit  equation. 

For  a  rotating  heat  pipe  with  the  physical  characteristics  as 
shown  in  Table  I,  the  amount  of  heat  that  can  be  transferred  from  the 
rotating  heat  pipe  is  limited  by  the  condensing  limit.  However,  the 
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limitations  imposed  by  the  sonic  limit,  boiling  limit,  and  entrainment  limit 
may  become  important  as  the  heat  pipe  geometry  and  operating  conditions 
vary. 

TABLE  I.  ROTATING  HEAT  PIPE  SPECIFICATIONS 


To  enhance  the  heat  transfer  capacity  of  the  rotating  heat  pipe, 
internally  finned  condensers  have  been  used  to  raise  the  condensing  limit 
line.  Thinner  films  occur  near  the  ridges  of  the  fins  while  thicker  films 
occur  in  the  troughs.  The  thinner  film  on  the  ridges  provides  a  lower 
thermal  resistance  to  heat  flow,  while  the  thicker  film  in  the  trough 
provides  a  higher  resistance.  A  compromise  between  the  improvement  on 
the  ridges  and  the  degradation  in  the  troughs  is  necessary  for  an  overall 
heat  transfer  improvement  [Ref.  5]. 

C.  ANALYSIS  OP  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 

Schafer  [Ref.  6]  developed  an  analytical  model  for  a  heat  pipe  with 
a  triangular  fin  profile  (figure  2).  This  model  was  developed  in  order  to 
raise  the  condensing  limit  by  the  addition  of  internal  fins.  An  assumption 
of  one-dimensional  heat  conduction  through  the  wadi  and  fin  was  made  for 
Schafer's  model. 

A  two-dimensionad  heat  conduction  model  using  a  Finite  Element  Method 
was  developed  for  this  same  case  by  Corley  [Ref.  7].  A  parabolic 
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temperature  distribution  along  the  fin  surface  was  assumed  by  Corley  [Ref. 
7].  A  significant  improvement  in  the  predicted  heat  transfer  performance 
was  indicated  by  his  res\alts.  By  using  the  two-dimensional  model  an 
increase  of  approximately  75  percent  in  the  heat  transfer  performance  was 
seen  over  the  results  from  the  use  of  the  one  dimensional  model.  However, 
Corley  [Ref.  7]  noted  that  at  the  fin  apex  an  error  as  great  as  50  percent 
was  possible  which  could  restdt  in  the  total  heat  transfer  being  in  error 
as  much  as  15  percent. 

A  modification  was  made  to  Corley’s  computer  program  by  Tantrakul 
[Ref.  8].  In  order  to  minimize  the  heat  transfer  error  at  the  apex  of  the 
fin  Tantrakul  increased  the  number  of  finite  elements  used.  His  results 
with  this  modification  converged  with  the  results  of  Corley. 

Purnomo  developed  a  linear  triangular  finite  element  model  (figure  3) 
used  in  a  two-dimensional  Finite  Element  Method  solution.  Purnomo's  [Ref. 
1]  Finite  Element  Method  program  also  worked  and  converged.  To  maximize 
the  heat  transfer  from  the  rotating  heat  pipe  the  condenser  geometry  was 
varied.  Using  Purnomo's  code  parametric  studies  were  conducted.  However, 
the  best  geometry  was  not  indicated  in  these  studies.  Purnomo's  code  was 
written  to  perform  one  analysis  at  a  time.  Davis  [Ref.  9]  modified 
Purnomo's  code  to  allow  for  numerous  aneilysis  to  be  made  using  the 
optimization  code  COPES/CONMIN.  Davis'  Finite  Element  .  Method  code 
incorporating  the  optimization  worked  and  converged,  resulting  in  an 
optimum  design  for  an  internally  finned  rotating  heat  pipe. 
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e  •  thleicaaas 
b  •  halghe  of  eh*  fia 


#  ■  coadaaaar  half  angl* 
a  ■  half  of  eh*  fia  vldeh 
c«  half  of  eh*  ecough  wtdeh 


Flg\ir«  3.  Condenser  Geometry  Considered  with  40  Linear  Triangul2u: 
Finite  Elements. 


D.  THESIS  OBJECTIVES 

The  objectives  of  this  thesis  were: 


1.  To  modify  Davis'  [Ref.  9]  computer  program  so  that  it  is  compatible 
with  the  ADS  (Automated  Design  Synthesis)  program  [Ref.  10]  and 
can  be  used  for  analysis  and  automated  design  of  rotating  heat 
pip«s. 
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2.  To  use  the  resulting  program  to  obtain  an  optimum  design  for  an 
internally  finned  rotating  heat  pipe  to  obtain  experimental  data  to 
compare  with  the  analytical  results. 

3.  To  use  the  resulting  program  to  obtain  numerical  results  in  place  of 
data  obtained  from  costly  experimental  operations. 


11 


n.  NUMERICAL  OPTIMIZATION 


A.  BACKGROUND 

The  parameter  that  ie  minimized  or  maximized  during  the  design 
process  is  called  the  design  objective.  The  design  objective  is  minimized 
or  meudmized  by  changing  the  design  variables  within  the  design  constraint 
limitations.  This  process  is  called  numerical  optimization.  An  assortment  of 
physical,  aesthetic,  economic  and,  on  occasion,  political  limitations  must  be 
met  by  the  design  constraints  for  the  design  to  be  acceptable.  For  the 
optimization  process  to  work,  the  design  criteria  must  be  described  in 
numerical  terms.  This  is  not  always  easy. 

A  computer  program  can  be  written  to  perform  tedious  and  repetitive 
calcvdations  necessary  to  optimize  the  problem  once  it  is  stated  in 
numerical  terms.  For  this  reason,  computer  analysis  is  commonplace  in  most 
engineering  organizations.  For  example,  in  heat  transfer  design  the 
configuration,  materials,  and  method  of  heat  removal  may  be  defined  and 
a  finite  element  analysis  computer  code  is  used  to  calculate  temperatures, 
heat  transfer  rates,  and  other  response  quantities  of  interest.  If  any  of 
these  parameters  are  not  within  prescribed  bounds,  the  engineer  may 
change  the  method  of  cooling  or  other  defined  quantity  and  rerun  thp 
program.  The  engineer  makes  the  actual  design  decisions,  the  computer 
code  only  provides  the  analysis  of  a  proposed  design.  This  is  the  commonly 
used  approach  which  is  called  computer-aided  design. 
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Analysis  codes  are  commonly  used  for  tradeoff  studies.  For  example, 
an  analysis  code  might  be  run  on  the  distance  a  truck  can  go  on  a  tank 
of  fuel.  For  different  loads,  different  distances  are  calculated  which  can 
be  used  in  a  range-payload  study. 

PviUy  automated  design  is  the  logical  next  step  to  computer-aided 
design.  The  computer  makes  the  actual  design  decisions  or  trade-off 
studies  based  on  input  criteria  in  fully  automated  design.  Minimal 
information  is  requested  from  the  operator  during  the  act.ual  design 
process.  Numerical  optimization  offers  numerous  improvements  over  the 
traditional  approach  to  design.  These  improvements  include:  time  reduction 
in  design  decision  making;  a  rational,  directed  design  procedure;  and  the 
procedure  is  unbiased  by  intuition  or  experience.  The  probability  of 
obtaining  a  non -traditional  solution  is  thereby  improved.  Engineering 
intxution  and  experience  are  still  necessary  to  decide  if  the  design  obtained 
is  an  improvement  and  feasible. 

B.  AUTOMATED  DESIGN  SYNTHESIS  (ADS) 

Vanderplaats  [Ref.  10]  developed  a  general  purpose  numerical 
optimization  program  containing  a  variety  of  algorithms,  ADS.  ADS  is  a 
FORTRAN  program  that  optimizes  a  numerically  defined  objective  function 
subject  to  a  set  of  constraint  limits.  The  solution  of  the  problem  is 
separated  into  three  levels: 

1.  Strategy  -  Optinization  strategy  such  as  Augmented  Lagrange 
Multiplier  method  or  Sequential  Linear  Programming. 

2.  Optimizer  -  Actual  algorithm  to  perform  the  optimization. 

3.  One-Dimensional  Search  -  Line  search  routine  used  by  optimizer. 
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Flexibility  to  solve  a  wide  variety  of  engineering  design  problems  is 
given  by  the  combinations  of  nine  strategies,  five  optimizers,  and  eight 
one-dimensional  search  options.  The  following  definitions  are  necessary  to 
discuss  the  use  of  ADS: 


1.  Design  Variables  -  Those  parameters  which  the  optimization  program 
is  permitted  to  change  within  allowed  bounds  in  order  to  improve 
the  design.  Design  variables  appear  only  on  the  right  hand  side  of 
an  equation  and  are  continuous. 

2.  Design  Constraints  -  An  inequality  constrzunt  reqrures  that  some 
function  of  the  design  variable(s)  remain  less  than  a  specified  value. 
Design  constraints  may  be  linear  or  nonlinear,  implicit  or  explicit, 
but  they  must  be  continuous  functions  of  the  design  variable. 

3.  Objective  Function  -  The  paurameter  which  is  going  to  be  minimized 
or  maximized  during  the  optimization  process.  The  objective  function 
may  be  linear  or  nonlinear,  implicit  or  explicit,  and  must  be  a 
continuous  function  of  the  design  variables.  The  objective  function 
usually  appear  on  the  left  side  of  an  equation. 


C.  PROGRAMMUIG  GUIDELINES 

Any  computer  code  developed  for  engineering  analysis  should  be 
written  in  such  a  way  that  it  is  easily  coupled  to  a  general  purpose 
optimization  program  such  as  ADS.  Therefore,  a  general  programming 
practice  is  outlined  here  which  in  no  way  inhibits  the  use  of  the  computer 
program  in  its  traditional  role  as  an  analytical  tool,  but  allows  for  simple 
adaption  to  ADS. 

ADS  is  called  by  a  user-supplied  calling  program.  ADS  does  not  rail 
any  user-supplied  subroutines.  Instead,  ADS  returns  control  to  the  caUing 
program  when  function  or  gradient  information  is  needed.  The  required 
information  is  evaluated  and  ADS  is  called  again.  This  provides  considerable 
flexibility  in  program  organization  and  restart  capabilLties.  Various  internal 


parameters  are  defined  on  the  first  call  to  ADS  which  work  well  for  the 
"average"  optimization  task.  However,  it  is  often  desirable  to  change  these 
in  order  to  gain  maximum  utility  of  the  program.  Figure  4  is  the  program 
flow  diagram  for  the  case  where  the  user  wishes  to  over-ride  one  or  more 
internal  parameters,  such  as  scaling,  convergence  criteria,  or  maximum 
number  of  iterations. 

After  initialization  of  basic  parameters  and  arrays,  the  information 
parameter,  INFO,  is  set  to  -2.  ADS  is  then  called  to  initialize  all  internal 
parameters  and  allocate  storage  space  for  internal  arrays.  Control  is  then 
returned  to  the  user,  at  which  point  these  parameters,  for  example 
convergence  criteria,  can  be  overridden  if  desired.  At  this  point,  INFO  will 
have  a  value  of  1  and  the  user  must  evaluate  the  objective  function,  OBJ, 
and  constraint  functions.  ADS  is  called  again  and  the  optimization  proceeds. 
Since,  in  this  case,  the  gradient  calculation  control,  IGRAD,  has  a  value  of 
zero,  all  gradient  information  is  calculated  by  finite  difference  within  ADS. 
When  INFO  has  a  value  of  zero,  optimization  is  complete. 
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BEGIN 


DIMENSION  ARRAYS 
DEFINE  BASIC  VARIABLES 
IGRAD=0  (USE  FINITE  DIFFERENCE  GRADIENTS) 

INFO  =  -2 
CALL  ADS  (INFO...) 

IF  INFO  =  0,  EXIT.  ERROR  WAS  DETECTED 
ELSE 

OVER-RIDE  DEFAULT  PARAMETERS  IN  ARRAYS  WK  AND  IWK  IF 

DESIRED 

CALL  ADS  (INFO...) 


NO 


YES 


INFO  =  0 


EVALUATE  OBJECTIVE  EXIT  OPTIMIZATION 

AND  CONSTRAINTS  IS  COMPLETE 


Figure  4.  Program  Flow  Logic:  Over-Ride  Default  Parameters,  Finite 
Difference  Gradients  10]* 
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m.  FINITE  ELEMENT  SOLUTION 

A.  REVIEW  OF  THE  PREVIOUS  ANALYSIS 

As  stated  previously,  the  heat  transfer  solution  for  a  one-dimensional 
model  of  an  internally  finned  rotating  heat  pipe  was  studied  by  Schafer 
[Ref.  6],  The  two-dimensional  model  was  studied  by  Corley  [Ref.  7].  The 
same  assumptions  and  boundary  conditions,  similar  to  those  used  in  the 
Nusslet  analysis  of  film  condensation  on  a  flat  wall,  and  based  upon  the 
analysis  of  Ballback  [Ref.  2]  were  used  for  both.  The  more  important  of 
these  assumptions  are: 

1.  steady  state  operation, 

2.  film  condensation,  as  opposed  to  dropwise  condensation, 

3.  laminar  flow  of  the  condensate  film  along  both  the  fin  and  the 
trough, 

4.  static  balance  of  forces  within  the  condensate, 

5.  one-dimensional  conduction  heat  transfer  through  the  film  thickness 
(no  convective  heat  transfer  in  the  condensate  film), 

6.  no  liquid -vapor  interfacial  shear  forces, 

7.  no  condensate  subcooling, 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of  the  wall  serrtion 
(symmetry  conditions),  as  shown  in  figure  5, 

9.  saturation  temperature  at  the  fin  apex, 

10.  zero  film  thickness  at  the  fin  apex,  and 

11.  negligible  curvature  of  the  condenser  wall. 
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Figure  3  shows  the  linear  triangular  finite  element  model  developed 
by  Purnomo  [Ref.  1]  for  use  in  obtaining  a  two-dimensional  Finite  Element 


solution. 

The  assumption  that  was  used  by  Corley  [Ref.  7]  that  the  fin  apex  was 
at  the  saturation  temperature  of  the  working  fluid  was  modified  by 
Purnomo  [Ref.  1].  The  value  of  the  temperature  at  the  apex  of  the  fin  was 
adlowed  to  float  and  a  parabolic  temperature  distribution  was  assumed  along 
the  fin  surface. 

Purnomo's  problem  statement  for  the  formulation  of  the  Finite  Element 
Method  as  shown  in  figure  5  is: 


dX^  dy‘ 


-0 


with  the  following  -boundary  ccndidons: 


(6) 


1.  along  boundary  1,  -K  dT]dn-h^a'-Tsat) 

2.  along  boundary  3,  -K  dTjdn-hJJ-T’o) 


3.  along  boundaries  2  and  4, 


dn 


-0 


A  detailed  description  of  the  numerical  formulation  is  presented  in  his 
thesis. 

Davis  used  Constrained  Function  Minimization  (CON  MIN)  as  an 
optimization  program.  CONMIN  is  a  FORTRAN  program  in  subrovitine  form. 


cona-tauit 


% 

iix  *  ill  .  0 

3x*  3y^ 

a)  -  :<  =  h,  (T  -  Along  Boundary  Cl] 

d  n  L  s«w 

b)  -  ^  s  h2  CT  -  T^)  Along  Boundary  C3] 

(-)  3  a  Along  Boundarias  C2]  and  C'*’] 


Elgure  5.  Differential  Equation  and  Boundaury  Conditions  Considered 
in  the  Analysis  of  Purnomo  [Ref.  1]. 


Vanderplaats  [Ref.  11]  developed  the  Control  Program  for  Engineering 
Synthesis  (COPES)  as  a  main  program  to  simplify  the  use  of  CONMIN.  Davis' 
computer  program  was  written  in  subroutine  form  with  SUBROUTINE  ANALIZ 
(ICALC)  2is  the  main  routine.  The  name  ANALIZ  is  compatible  with  the  COPES 
program  and  ICALC  is  a  calculation  control.  Subroutine  ANALIZ  calls  other 
subroutines  as  needed: 

1.  the  routine  "COORD"  used  to  define  positions  of  system  coordinate 
points, 

2.  the  routine  "FORMAF"  used  to  formulate  the  Finite  Element  Method 
equations, 

3.  the  routine  "BANDEC"  as  an  equation  solver  for  a  symmetric  matrix 
which  has  been  transformed  into  banded  form,  and 

4.  the  routine  "DPLORT"  used  to  compute  the  roots  of  a  real  polynomial 
using  a  Mewton-Raphson  derivative  technique. 

B.  THE  COMPUTER  PROGRAM  DEVELOPMENT 

The  basis  for  the  present  analysis  code  is  Davis'  [Ref.  9]  two- 
dimensional  finite  element  program.  Davis'  code  was  checked  for  validity  as 
the  first  task  undertaken  in  the  development  of  this  thesis.  An  error  was 
discovered  in  calcxilating  the  fin  condensate  film  thickness  (AZS).  In  the 
initial  calculation  of  HDEN  only,  the  cubed  term  was  merely  multiplied  by 
three.  The  correct  form  of  the  equation  is  shown  below: 

HDEN  - 

The  effect  of  this  error  was  minimal  since  subsequent  equations  were 
correct,  a  0.00016%  difference  in  the  condensate  level  was  noted. 
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The  next  task  undertaken  was  to  adapt  the  analysis  code  to  permit 
automated  design  analysis  using  ADS.  Many  modifications  were  made,  some 
of  which  are  mentioned  here.  The  program  was  rewritten  to  include  a  main 
program  from  which  ADS  is  called  and  subroutines  to  perform  various 
mathematical  functions.  Subroutine  ANALIZ  was  deleted  since  the  double 
precision  version  of  ADS  (DADS)  was  used.  Modifications  were  also  made  to 
generalize  the  code  and  minimize  the  changes  needed  when  varying  the 
number  of  finite  elements  used. 

A  listing  of  the  revised  computer  program  is  included  as  the 
Appendix. 

C.  DESIGII  OPTIMIZATIOR 

There  are  thirteen  parameters  that  can  be  used  as  design  variables. 
There  are  geometric  or  functional  parameters  of  the  rotating  heat  pipe  or 
the  properties  of  the  working  fluid  or  environment.  The  possible  design 
variables,  possible  constr2dnt  functions,  and  the  objective  function  are 
listed  below  in  Fortran.  This  code  can  pursue  a  wide  variety  of  design 
problems. 

The  addition  of  fins  by  the  designer  increases  the  surface  area  which 
increases  the  heat  transfer  rate  through  the  condenser  wall.  However,  the 
addition  of  fins  decreases  the  cross-sectional  area  through  each  fin  for 
conduction  and  decreases  the  trough  width  which  increases  the  condensatp 
thickness  in  the  trough.  The  increased  condensate  thickness  decreases  the 
heat  transfer  and  if  increased  to  the  point  of  covering  the  fins  it  could 
dramatically  reduce  the  heat  transfer  through  the  fin. 


21 


TABLE  n.  DESIGN  VARIABLES 


DESIGN  VARIABLES _ 

BFIN  (fin  height) 

CANGL  (cone  half  angle) 

CLI  (condenser  length 
PANGL  (fin  half  angle) 

HINF  (external  convective  heat  transfer  coefficient) 
NPIN  (number  of  fins) 

R2I  (intermediate  radius) 

RBASEI  (inside  radius  of  condenser  base) 

RMP  (rotational  speed  of  the  heat  pipe) 

THICKI  (condenser  waU  thickness) 

TINF  (external  temperature) 

TSS  (saturation  temperature) 

TZ  (nodal  point  temperature) _ 

CONSTRAINT  FDNCTIONS _ 

BOA  (ratio  of  fin  height  over  fin  width) 

ZOA  (ratio  of  trough  width  over  fin  width) 

DEL(NI)  (condensate  thickness) _ 

OBJECTIVE  FUNCTION 


OBJ 


1 

4 


,i 

i 

(I 


li 


< 

•i 


■I 


The  purpose  of  the  design  study  was  to  determine  the  fin  height, 
number  of  fins,  and  fin  half  angle  which  would  maximize  the  heat  transfer 
rate.  It  was  then  decided  that  the  design  variables  would  be  BFIN,  FANGL, 
ouid  NFIN.  The  number  of  fins  was  chosen  vice  the  ratio  of  trough  width 
to  fin  width  (ZOA)  as  was  previously  used.  This  decision  was  made  since 
it  is  easier  to  think  in  terms  of  the  number  of  fins  vice  a  ratio.  Other 
potential  design  variables  were  held  constant.  The  objective  function  to  be 
maximized  was  OBJ=QT+QTF,  the  heat  transfer  through  the  fin  plus  the  heat 
transfer  through  the  trough. 

The  code  was  run  with  the  three  design  variables  using  the  internal 
scalar  default  parameters  in  ADS.  The  objective  function  was  calcvilated 
using  the  input  values  of  the  design  variables.  The  ADS  program  then 
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changed  the  first  design  variable  keeping  the  other  two  design  variables 
at  the  input  values.  The  calculations  were  made  aging  yielding  a  different 
value  for  the  objective  function.  The  new  value  would  then  be  compared 
to  the  previous  value  of  the  objective  function.  If  the  difference  in  the 
objective  function  was  not  greater  then  the  internal  scalar  default  value, 
the  first  design  variable  would  be  returned  to  its  original  value  and  the 
process  repeated  with  the  second  design  variable.  If  the  difference  in  the 
objective  function  was  still  not  greater  then  the  internal  scalar  defavilt 
value,  the  second  design  variable  would  be  returned  to  its  initial  value  and 
the  process  repeated  with  the  third  design  variable. 

Each  time  the  program  was  run,  the  optimization  code  would  choose 
BFIN  as  the  design  variable  to  change  first  as  it  had  the  greater  influence 
on  the  objective  function.  The  remaining  two  variables  would  then  either 
be  kept  constant  or  the  number  of  fins  would  be  maximized  and  the  fin 
half  angle  minimized.  Consistent  results  were  not  obtained  with  this  method. 

To  improve  the  results,  the  internal  scalar  parameters  were  modified. 
These  modifications  included  the  constraint  tolerances,  the  absolute  and 
relative  convergence  criteria,  the  absolute  and  relative  change  in  the 
design  variables,  the  absolute  and  relative  change  in  the  objective 
function,  the  minimum  absolute  value  of  the  finite  difference  step  when 
calculation  gradients,  and  the  initial  relative  move  limit.  Better  results  were 
obtained  as  seen  in  the  higher  value  for  the  objfjctive  function.  However, 
the  results  were  still  not  consistent  depending  on  the  initial  values  used 
for  the  three  design  variables.  At  this  point,  it  was  decided  to  concentrate 
on  one  design  variable  and  on  the  basis  of  the  previous  calculations  thp 
design  variable  chosen  was  BFIN. 
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The  external  surface  temperature  was  set  equal  to  the  working  fluid 
saturation  temperature  and  the  theoretical  upper  limit  on  the  heat  transfer 
was  calciilated  for  comparison.  An  assumption  is  made  that  there  is  no 
thermal  resistance  across  the  condensate  or  the  condenser  wall.  The  upper 
limit  of  the  heat  transfer  rate  was  predicted  to  be  69,492  BTU/HR  using  the 
following  formula: 

o) 

where 

h  =  outside  convective  heat  transfer  coefficient  (5000  BTU/HR ’FT^'F) 
r  =  average  outside  radius  of  condenser  wall  (0.07373  FT) 

1  =  condenser  length  (0.75  FT) 

Tyj’i  s  temperature  of  the  outside  wall  (100*F) 

e 

T*  a  ambient  temperature  (60*P) 

Based  on  engineering  judgement  certain  constraints  were  placed  on  the 
design.  These  constraints  were  applied  to  the  number  of  fins  (not  to 
exceed  400)  and  the  minimum  fin  half  angle  (not  to  be  less  than  10 
degrees).  The  values  were  based  on  structural  and  manufacturing 
considerations. 
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IV.  RESULTS 


A.  INTRODUCTION 

The  purpose  of  the  design  optimization  was  to-  maximize  the  heat 
transfer  rate.  This  was  accomplished  by  using  the  computer  code  in 
conjunction  with  ADS.  Numerical  results  are  discussed  below. 

B.  CONSTRAINED  OPTIMIZATION 

In  the  design  problem  undertaken  to  determine  the  optimum  internal 
geometry  for  the  maximum  heat  transfer,  numerous  runs  were  made  for  a 
condenser  made  of  copper.  This  material  has  a  thermal  conductivity  of  230 
BTU/HR’FT’P.  The  working  fluid  was  water. 

Since  the  fin  height  (BFIN)  was  the  design  variable,  the  initial  runs 
investigated  whether  there  vas  an  optimum  fin  height.  Initially  the  fin  half 
angle  was  held  constant  at  10  degrees  and  the  number  of  fins  was  varied 
from  150  to  400.  In  each  case,  for  the  optimum  design,  the  fin  height  was 
maximized  and  the  trough  was  eliminated,  as  seen  in  figures  6  and  7. 

The  number  of  fins  was  then  held  constant  and  the  fin  half  angle  was 
varied  from  10  to  25  degrees  with  the  fin  height  remaining  the  design 
variable.  Once  again,  the  greatest  heat  transfer  rate  was  achieved  with  the 
highest  fin  height  for  each  number  of  fins  (figure  8., 

As  seen  in  figure  9,  the  highest  heat  transfer  rate  achieved  was  for 
400  fins  with  a  10  degree  fin  half  angle  and  a  fin  height  of  0.0345  inches. 
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Figure  6.  Heat  Transfer  Rate  vs.  Fin  Height. 
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Figure  7.  Heat  Transfer  Rate  vs.  Pin  Height. 
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Fiqxire  8.  Heat  Transfer  Rate  vs.  Pin  Half  Angle  (Optimum  Fin  Height). 


Figure  10  shows  a  plot  of  heat  transfer  rate  versus  fin  hadf  angle  for 
a  condenser  with  between  100  and  400  fins.  The  heat  transfer  rate,  as  a 
function  of  fin  half  angle  for  a  constant  fin  height,  increases  as  the  fin 
half  angle  increases.  Davis  [Ref.  9]  concluded  that  the  heat  transfer  rate 
increased  with  an  increase  in  fin  half  angle.  This  is  correct  if  the  fin 
height  is  kept  constant.  As  the  fin  half  angle  increases,  the  surface  area 
of  the  fin  increases.  The  added  surface  area  also  has  a  thinner  film  of 
condensate  on  it  which  offers  lower  thermal  resistance.  The  trough  area 
decreases  and  the  condensate  film  in  the  trough  thickens,  increasing  the 
thermal  resistance.  This  degradation  does  not  offset  the  gain  in  the  heat 
transfer  rate  caused  by  the  fin. 

For  external  heat  transfer  coefficients  varying  from  1000-50,000 
BTU/HR«FT^«P,  the  same  optimum  design  geometry  for  a  maximum  heat 
transfer  rate  was  obtained,  which  is  stated  in  table  III  below.  Figure  11 
shows  the  strong  influence  the  external  heat  transfer  coefficient  has  on 
the  heat  transfer  rate. 

In  figure  12,  the  effect  of  the  rotating  speed  on  the  heat  transfer 
rate  is  seen.  As  the  rotational  speed  increases,  the  heat  transfer  rate 
increases,  this  is  caused  by  an  increase  in  the  element  heat  transfer 
coefficient  and  a  decrease  in  the  condensate  thickness  on  the  fin  which 
lowers  the  thermal  resistance. 
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Figure  11.  Heat  Transfer  Rate  vs.  Heat  Transfer  Coefficient. 
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Figxire  13.  Condensate  Level  vs.  Position  (100-400  Fins). 
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When  the  number  of  fins  was  increased  from  100  to  400,  maintaining 
the  same  fin  half  angle,  the  condensate  level  decreased  with  the  increased 
heat  transfer  rate  (figure  13).  This  occurred  because  the  thinner  film  over 
the  fins  decreased  the  resistance  across  the  film  which  raised  the 
temperatures  along  the  fin,  which  in  combination  with  the  lower  height  fins 
increased  the  temperature  on  the  outside  of  the  pipe.  This  increase  in 
temperature  brings  the  operation  closer  to  the  condensing  limit.  When  the 
fin  half  angle  was  increased  for  r.  specified  number  and  height  of  fins,  the 
condensate  level  decreased  due  to  the  increeised  trough  width  (figure  14). 


TABLE  m.  OPTIMIZATION  RESULTS  FOR 
VARIOnS  HEAT  TRANSFER  COEFFICIENTS 


OPTIMUM 

1  DESIGN 

Fin  Height 

0.0345  inches 

Fin  Half  Angle 

10.0  degrees 

Number  of  Fins 

400 

^external 

HEAT  TRANSFER  RATE 

(BTU/HR«FT^«F) 

(BTU/HR) 

1000 

13765 

5000 

62252 

50000 

261003 

Figures  15  and  16  show  the  effect  the  increase  of  surface  area  has 
on  the  heat  transfer  rate.  In  figure  15,  the  fin  height  for  400  fins  with  a 
10  degree  half  angle  is  plotted  against  the  heat  transfer  rate.  As  the  fin 
height  is  increased  up  to  a  maximum  value  of  0.0345  inches  the  resulting 
design  is  a  sawtooth.  Figure  16  shows  the  effect  adding  more  fins  has  on 
the  heat  transfer  rate  for  a  constant  height  fin.  The  greatest  increase  in 
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the  heat  transfer  rate  is  seen  from  the  addition  of  100  fins  from  a  smooth 


tube. 

In  figure  17,  the  ratio  of  the  actual  surface  area  over  the  surface 
aurea  for  a  smooth  tube  is  plotted  versus  the  number  of  fins.  The  fin 
height  and  the  fin  half  angle  are  both  held  constant.  As  expected,  the 
surface  area  ratio  increases  in  a  linear  manner  as  the  number  of  fins 
increases.  The  ratio  of  the  heat  transfer  rate  over  the  heat  transfer  rate 
for  a  smooth  tube  is  plotted  for  two  different  heat  transfer  coefficients, 
h=1000  BTU/HR«FT^»F  and  5000  BTU/HR»FT^»F.  An  increase  in  the  number 
of  fins  results  in  not  only  aui  increase  in  the  area  ratio  but  also  an 
increase  in  the  heat  transfer  ratio.  The  increase  in  the  heat  transfer  ratio 
is  greatest  when  going  from  a  smooth  tube  to  a  tube  with  100  fins.  The 
heat  transfer  ratio  increase  is  greater  for  the  heat  transfer  coefficient 
equal  to  5000  BTU/HR«FT^«F.  This  is  because  the  heat  transfer  coefficient 
has  a  direct  effect  on  the  heat  transfer  rate,  that  is, 

Q  -  M(T-T^ 

Both  curves  approach  an  asymptotic  value.  However,  the  curve  with  the 
lower  heat  transfer  coefficient  approaches  this  aisymptotic  value  with  a 
fewer  number  of  fins.  Additionally,  in  view  of  the  relatively  small  increase 
in  the  heat  transfer  ratio  by  the  addition  of  fins  for  the  lower  heat 
transfer  coefficient,  consideration  should  be  given  to  the  cost  of 
manufacturing  the  fins  versus  the  benefit  derived  by  their  addition. 
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Figure  14.  Condensate  Level  vs.  Position  (400  Pins,  10-25  Degree  Fin 
Half  Angle). 


NUMBER  OF  RNS 


Pigiiure  16.  Heat  Transfer  Rate  vs.  Number  of  Pins  (Constant  Pin 
Height). 
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figure  17.  Heat  Transfer  and  Area  Ratios  vs.  Number  of  Pins. 
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C.  AUTOMATED  DESIGN  SYNTHESIS  (ADS) 

This  optindzatiDn  project  was  done  on  the  IBM  mainframe  using  the 
ADS  optimization  program.  As  stated  previously,  ADS  is  a  general  purpose 
numerical  optimization  program  with  a  variety  of  algorithms  that  can  be 
used  to  tailor  the  solution.  The  solution  is  separated  into  three  levels  in 
ADS:  strategy,  optimizer,  and  one-dimensional  search.  For  this  problem  the 
following  combination  of  algorithms  were  used: 

1.  Strategy:  Sequential  Linear  Programming 

2.  Optimizer:  Modified  Method  of  Feasible  Directions  for  constrained 
minimization 

3.  One-Dimensional  Search:  Golden  Section  Method  followed  by 
polynomial  interpolation 

The  strategy  used  linearizes  a  nonlinear  problem  by  a  first  order 
Taylor  series  expansion  of  the  objective  and  constraint  functions.  The 
solution  to  this  linear  approximation  is  obtained.  The  problem  is  linearized 
again  about  this  point  and  the  new  problem  is  solved  with  the  process 
being  repeated  until  a  precise  solution  is  achieved. 

The  optimizer  chosen  is  used  to  find  a  search  direction  which  will 
minimize  the  objective  function  while  maintaining  feasibility. 

The  combination  of  strategy,  optimizer,  and  one-dimensional  search 
chosen  is  not  the  only  one  available,  nor  is  it  necessarily  the  most 
efficient.  It  did  yield  results  that  were  maximized  and  were  within  the 
constraint  tolerances. 

The  ADS  optimization  program  is  complicated  by  the  numerous  internal 
parameters  which  must  be  changed  to  obtain  an  optimal  design. 
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Complications  aurise  when  there  is  a  vast  difference  in  the  scales  of  the 
design  variables.  The  design  variables  themselves  must  be  scaled  which  is 
further  complicated  when  the  variable  itself  covers  a  wide  range.  ADS  also, 
in  certain  cases,  allows  for  constraints  to  be  violated.  In  some  instances 
this  might  be  acceptable  but  not  in  this  case. 

AOS  does  not  have  a  scoping  mechanism,  that  is  the  ability  to 
decrease  the  rate  of  change  of  the  design  variable,  and  therefore  to  obtain 
a  precise  answer  the  internal  paurameters  must  be  changed  repeatedly.  ADS 
also  does  not  recognize  integers,  all  numbers  are  real  therefore,  depending 
on  the  answer  given,  it  may  be  necessary  to  round  up  or  down. 
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V.  CONCLUSIOIfS 


1.  For  an  independent  increase  in  fin  half  angle,  rotational  speed  or 
number  of  fins  an  increase  is  seen  in  the  heat  transfer  rate.  As  the 
parameters  increase,  the  heat  transfer  rate  levels  off  at  the  theoretical 
maximum  heat  transfer  rate  for  the  heat  pipe.  A  decrease  in  the  fin  half 
angle  with  a  corresponding  increase  in  the  fin  height  increases  the  heat 
transfer  rate.  If  the  fin  half  angle  is  decreased  while  the  fin  height  is 
kept  constant,  then  the  heat  transfer  rate  decreases. 

2.  Maximum  heat  transfer  occtirs  for  the  same  fin  geometry 
regardless  of  the  external  heat  transfer  coefficient.  For  a  specific 
condenser  radius,  as  many  fins  as  possible  should  be  machined  with  a 
minimum  fin  half  angle  at  the  maximum  fin  height. 

3.  The  computer  code  cam  be  used  for  single  analysis  or  the 
automated  design  of  an  internally  finned  rotating  heat  pipe. 

4.  The  benefit  of  adding  fins  is  dependent  on  the  external  heat 
tramsfer  coefficient.  Consideration  of  the  cost  of  mamufacturing  the  fins 
versus  the  increase  in  the  heat  transfer  rate  should  be  made. 
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VI.  RECOMMENDATIONS 


1.  Analyze  different  shaped  fins  including  rectangular  and  curved. 

2.  Modify  the  code  to  allow  for  silmutaneous  variations  of  more  than 
one  variable. 

3.  Use  different  working  flxuds  and  heat  pipe  materials  to  see  if  a 
different  internal  geometry  occurs  for  the  maiximum  heat  transfer  rate. 

4.  Modify  the  code  to  use  the  DOT  optimization  program  vice  ADS. 
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APPENDIX:  PROGRAM  LISTING 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


PROGRAM  OLSON 

It  * 

*  ANALYSIS  OF  ROTATING  HEAT  PIPE  ,  USING  TRIANGULAR  * 

*  ELEMENT  MODEL  * 

*  COMPILED  BY  MAJOR  IGNATIUS . S . PURNOMO  IN  JUNE  1978  * 

*  * 

*  MODIFIED  TO  PERMIT  NUMERICAL  OPTIMIZATION  * 

*  USING  COPES/CONMIN  * 

*  BY  LCDR  WILLIAM  A.  DAVIS,  JR.  IN  SEPTEMBER  1980  * 

*  MODIFIED  TO  PERMIT  USE  OF  THE  ADS  OPTIMIZATION  * 

*  ROUTINE  BY  LT.  G.L.  OLSON  IN  JUNE  1992  * 

* 

*  * 


CHARACTER* 20  NAME 

* 

COMMON/ADS/DF(21) ,G(10) ,1DG(100) ,IGRAD,INFO,IOPT,IONED,IPRINT, 
!ISTRAT,IWK(2000) ,IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W( 21 , 30 ) , WK( 5000 ) , 
; NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

* 

C 

COMMON/OLLIE/A( 200,50) ,AMTOT( 200) ,APS,B( 3 ) , BFIN , BOA, BVIN , C ( 3 ) , 
:CANGL,CF( 200) ,CK,CLI ,COF( 5) ,CW( 200) , DEL ( 200 ) ,DMDOT( 200),EA(3,3), 

; EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ( 100 ) ,H(200) , HINF , HZ ( 200 ) , 
:QB(200) ,QINC(200) ,QTINC(200) , QTOT , QTOTAL (1 0 0 ) ,R(200),RB(200), 
:RBASEI,R2I,RHOF(200) ,ROOTI(4) ,RPM, ROOTR( 4 ) ,T(200) ,TALFA,TB( 200 ) , 
:TCC(200) ,TE(200) , THICK, THICKI , TIB { 200 ) , TINF , TS ( 200 ) ,TSAT,TSS, 
:TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,ZOA, 
;DOBF,DOTH, ICOR( 200, 3) , IFF, JINT, JLC, JTC, KFF( 50) ,KFIN( 50 ) ,KT,NBAN, 
:NEL,NFIN,NSNP 

* 

* 

C  GUIDE  TO  FORTRAN  VARIABLE  NAMES 

C 


C 

AFOVAS 

FIN  AREA/SMOOTH  AREA 

c 

ALFA 

FIN  HALF  ANGLE  (RADIANS) 

c 

BFIN 

HEIGHT  OF  FIN  (INCHES) 

c 

BOA 

TANGENT  OF  THE  FIN  HALF  ANGLE 

c 

BVIN 

HEIGHT  OF  FIN  (FEET) 

c 

CAL  FA 

COSINE  OF  ALFA 

c 

CANGL 

CONE  HALF  ANGLE  (DEGREES) 

c 

CBASE 

INSIDE  CIRCUMFERENCE  OF  CONDENSER 

(  FEET 

) 

c 

CEXIT 

INSIDE  CIRCUMFERENCE  AT  CONDENSER 

EXIT 

(  FEET) 

c 

CF 

THERMAL  CONDUCTIVITY  OF  CONDENSATE 

FILM 

( BTU/HR  FT  F 

c 

CL 

CONDENSER  LENGTH  (FEET) 

c 

CLI 

CONDENSER  LENGTH  (INCHES) 

c 

CPHI 

COSINE  OF  PHI 

c 

CRIT 

CONVERGENCE  CRITERION 
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C  DEL  FILM  THICKNESS 

C  DF  G^IADIENT  OF  OBJECTIVE 

C  DIV  FLOATING  POINT  VALUE  OF  NDIV 

C  DMTOT  CONDENSATE  MASS  FLOW  RATE 

C  EPS  TROUGH  WIDTH  INCLUDING  INCREMENTAL  CHANGE 

C  EPSEX  TROUGH  WIDTH  AT  CONDENSER  EXIT 

C  EPSO  TROUGH  WIDTH  AT  START  OF  CONDENSER 

C  E2ERO  FIN  BASE  WIDTH 

C  F  FORCE  VECTOR  OF  SYSTEM 

C  FANGL  FIN  HALF  ANGLE  (DEGREES) 

C  G  CONSTRAINT  VALUES  ASSOCIATED  WITH  CURRENT  DESIGN 

C  H  CONVECTIVE  HEAT  TRANSFER  COEFFICIENT  ( BTU/HR  FT2  F) 

C  HFG  LATENT  HEAT  OF  VAPORIZATION  (BTU/LBM) 

C  IDG  CONSTRAINT  TYPE  IDENTIFIER 

C  I EL  THE  ELEMENT  NUMBER 

C  IFF  NO.  OF  ROWS  MINUS  ONE  OF  THE  UPPER  TRIANGULAR  FIN 

C  I FIN  EQUALS  0  FOR  COPPER,  AND  1  FOR  STAINLESS  STEEL 

C  I FLUID  EQUALS  0  FOR  WATER,  AND  1  FOR  FREON 

C  I GRAD  GRADIENT  CALCULATION  IDENTIFIER 

C  INFO  CONTROL  PARAMETER 

C  lONED  ONE  DIMENSIONAL  SEARCH  IDENTIFIER 

C  lOPT  OPTIMIZER  IDENTIFIER 

C  IPRINT  A  FOUR  DIGIT  PRINT  CONTROL 

C  ISTRAT  OPTIMIZATION  STRATEGY  IDENTIFIER 

C  JINT  NO.  OF  COLUMNS  PLUS  ONE  BELOW  TRIANGULAR  FIN 

C  JLC  NUMBER  OF  SYSTEM  NODAL  POINT  LOCATED  AT 

C  THE  CENTER  OP  SYSTEM  COORDINATE 

C  JTC  NUMBER  OF  SYSTEM  NODAL  POINT  LOCATED  AT 

C  THE  JUNCTION  OF  THE  SYMMETRY  BOUNDARY  AND 

C  THE  LINE  OP  INTERSECTION  BETWEEN  THE  FIN 

C  AND  THE  CONDENSER  WALL 

C  KPF  NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  ALONG 

C  THE  FIN  CONVECTIVE  BOUNDARY 

C  KFIN  NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  ON  THE 

C  SYSMMETRIC  BOUNDARY  OF  TRIANGULAR  FIN  SECTION 

C  NOT  COUNTING  POINTS  AT  BASE  AND  APEX 

C  KT  NUMBER  OF  COLUMNS  WITHIN  THE  TROUGH  WALL  SECTION 

C  NBAN  SYSTEM  BAND  WIDTH 

C  NBOTF  LAST  ELEMENT  AT  BOTTOM  SIDE 

C  NBOTI  FIRST  ELEMENT  AT  BOTTOM  SIDE 

C  NCON  NUMBER  OF  CONSTRAINTS 

C  NDOBF  NUMBER  OF  ROWS  WITHIN  THE  FIN 

C  NDOTH  NUMBER  OF  ROWS  WITHIN  THE  TROUGH 

C  NDIV  NUMBER  OF  INCREMENT 

C  NDV  NUMBER  OF  DESIGN  VARIABLES 

C  NEFB  ELEMENT  NUMBER  AT  BASE  OF  FIN 

C  NEL  NUMBER  OF  ELEMENTS 

C  NEST  ELEMENT  NUMBER  AT  END  OF  TROUGH 

C  NRA  NUMBER  OF  ROWS  IN  ARRAY  A 

C  NRWK  DIMENSIONAL  SIZE  OF  WK 

C  NSNP  NUMBER  OF  SYSTEM  NODAL  POINTS 

C  OBJ  VALUE  OF  THE  OBJECTIVE  FUNCTION  ASSOCIATED  WITH  X 

C  OMEGA  ROTATIONAL  SPEED  OF  HEAT  PIPE  (RAD/HR) 

C  PHI  CONE  HALF  ANGLE  (RADIANS) 

C  PI  PI 

C  R2I  DISTANCE  FROM  CENTERLINE  OF  THE  HEAT  PIPE  TO  HALF  THE  FIN 

C  HEIGHT 

C  RBASE  INSIDE  RADIUS  OF  CONDENSER  BASE  (FEET) 

C  RBASEI  INSIDE  RADIUS  OF  CONDENSER  BASE  (INCHES) 

C  REXIT  INSIDE  RADIUS  OF  CONDENSER  EXIT  (FEET) 
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»n*  »n*  *o**n»******oo 


c 

RPM 

REVOLUTIONS  PER  MINUTE 

c 

S 

VECTOR  OF  DESIGN  VARIABLES 

c 

SAL  FA 

SINE  OF  ALFA 

c 

SPHI 

SINE  OF  PHI 

c 

SURFAR 

SURFACE  AREA 

c 

THICK 

CONDENSER  WALL  THICKNESS  (FEET) 

c 

THICKI 

CONDENSER  WALL  THICKNESS  (INCHES) 

c 

TPHI 

TANGENT  OF  PHI 

c 

TZ 

AMBIENT  TEMPERATURE 

c 

UF 

VISCOSITY 

c 

VLB 

LOWER  BOUNDS  ON  THE  DESIGN  VARIABLES  ' 

c 

VUB 

UPPER  BOUNDS  ON  THE  DESIGN  VARIABLES 

c 

WK 

REAL  WORK  ARRAY 

c 

ZFIN 

NUMBER  OF  FINS 

c 

ZOA 

RATIO  OF  TROUGH  WIDTH  TO  FIN  BASE  WIDTH 

c 

ZSTAR 

SURFACE  LENGTH  OF  THE  FIN  MINUS  THE  SURFACE  LENGTH 

c 

COVERED  BY  THE  CONDENSATE  IN  THE  TROUGH 

c 

ZZERO 

SURFACE  LENGTH  OP  FIN 

C 

C 

PRINT*,  'INPUT  PILE  NAME' 

READ*,  NAME 

OPEN ( 1 0 , FILE-NAME ) 

OPEN( 15,FILE-'/HTPIPE  OUTPUT') 

OPEN (14, PILE- '/DUMP  OUTPUT') 

OPEN( 13, FILE- '/GRAPH  OUTPUT') 

* 

* 

THE  FOLLOWING  READS  INPUT  DATA,  PERFORMS  HEAT 
TRANSFER  ANALYSIS,  AND  PRINTS  RESULTS. 


*****  INPUT  MODE  ***** 


ELEMENT  CONNECTIVITIES 

READ  (10,420)  NEL,NSNP,NBAN, IFLUID, IFIN 
WRITE  (15,430)  NEL , NSNP , NBAN 
WRITE  (15,435)  IFLUID, IFIN 
WRITE  (15,436) 

WRITE. (15,437) 

READ  (10,440)  ( ICL , ( ICOR ( lEL , I ) , I-l , 3 ) , I EL-1 , NEL ) 
WRITE  (15,450) 

THE  CONDENSER  GEOMETRY 

READ  (10,460)  CLI , CANGL , RBASEI , R2I , THICKI , BFIN, TZ 
WRITE  (15,470)  CLI , CANGL , RBASEI , R2l , THICKI , BFIN , TZ 
READ  (10,480)  NDIV , NEST , NEFB , NBOTI , NBOTF 
WRITE  (15,490)  NDIV, NEST, NEFB, NBOTI , NBOTF 

DATA  FOR  RUNNING 

READ  (10,500)  RPM,TSS,TINF,HINF 
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* 

c 

* 


WRITE  (15,510)  RPM,TSS,TINF,HINF 
THE  CONVERGENCE  CRITERIAN 


READ  (10,520)  CRIT 
WRITE  (15,530)  CRIT 

* 

C  INTERNAL  FIN  GEOMETRY 

* 

READ  (10,540)  FANGL,NFIN 
WRITE  (15,550)  FANGL.NFIN 
WRITE(*,*)  FANGL,NFIN 
READ  (10,560)  IFF 
WRITE  (15,570)  IFF 

READ  (10,580)  ( KF1N( I ) , KFF( I ) , I-l , IFF ) 

READ  (10,590)  NDOBP,NDOTH, JTC, JLC, JINT,KT 

NHB-NEFB/2 

NBF-NBOTF+1 

DOBF  -  FLOAT (NDOBP) 

DOTH  -  FLOAT (NDOTH) 

WRITE  (15,600)  ICOR(NBOTI,2) ,ICOR(NEFB,l) ,ICOR(NEST,l) , 
lICOR(NBOTF,l) 

* 

*  SET  CONSTRAINTS 
NRA-21 
NCOLA-30 
NRWK-  5000 
NRIWK-2000 
NDV-  1 

NCON-  2 
I GRAD-  0 

* 

*  INITIAL  DESIGN 
S(l)-  BPIN 

* 

*  BOUNDS 

VLB(l)-  0.0000001 
VUB(l)-  0.75 

* 

C  IDENTIFY  CONSTRAINTS 

C  NONLINEAR  CONSTRAINT 

IDG(1)-1 

C  LINEAR  CONSTRAINT 

IDG(2)-2 

* 

PRINT*,  'INPUT  THE  VALUES  FOR  ISTRAT , lOPT , lONED  AND  IPRINT' 
READ*,  ISTRAT, lOPT, lONED, IPRINT 


C  INITIALIZE  COUNTER 

NO-0 . 0 

C  CHANGE  THE  INTERNAL  PARAMETERS 

INFO— 2 

CALL  DADS( INFO, ISTRAT, lOPT , lONED , IPRINT , IGRAD , NDV, NCON , S , VLB , 
:VUB,OBJ,G, IDG,NGT, IZ , DF , W, NRA , NCOLA, WK , NRWK , IWK , NRIWK ) 

IWK( 2 )-0 
IWK( 3 )-200 
IWK( 5)-4 
IWK( 7 )»500 
WK(3)— 0.5 
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nno*'  oo**n 


WK(6)— 0.01 

WK(8)-0.05 

WK(9)-0.50 

WK(10)-0.05 

WK(ll)-0.005 

WK(13)-0.001 

WK(14)-0.0001 

WK(21)-0.004 

WK(22)>0.002 

WK(26)«0.004 

WK(37)-0. 0000000001 

10  CALL  DADS(INFO,ISTRAT,IOPT,IONED,IPRINT,IGRAD,NDV,NCON,S,VLB, 
; VUB , OB J , G , I DG , NGT , I Z , D  F , W , NRA , NCOLA , WK , NRWK , I WK , NRI WK ) 

IF  (INFO.EQ.O)  GO  TO  360 

*****  EXECUTION  MODE  ***** 

NO-NO+1 

CONVERT  UNITS  OF  ALL  DIMENSIONAL  PARAMETERS 
TO  FEET.  CONVERT  UNITS  OF  ANGLES  TO  RADIANS. 

R2I-RBASEI- ( S ( 1 )/2 ) 

CL-CLI/12.0 
R2-R2I/12.0 
RBASE-RBASEI/1 2 . 0 
BVIN-S(1)/12.0 
DIV-FLOAT(MDIV) 

PI-3.14159265358979 
PHl-2 . 0*CANGL*PI/360 . 0 
SPHI-SIN(PHI) 

CPHI-COS(PHl) 

TPHI-TAN(PHI) 

DELX-CL/DIV 
CBASE-2.0*PI*RBASE 
REXIT-RBASE+CL*TPHI 
CEXIT-2 . 0’*PI*REXIT 
THICK-THICKI/12 . 0 
ALFA-FANGL*2 . 0*PI/360 . 0 
SALFA-SIN(ALFA) 

CALFA-COS(ALFA) 

TALFA-TAN(ALFA) 

EZERO-2 .0*(S(1)/12.0) *TALFA 
ZOA-( (CBASE-( EZERO*NFIN) )/NFIN)/EZERO 

BOUNDARY  CONDITIONS  AND  TEMPERATURE  ESTIMATES 
ALONG  THE  FIN  BOUNDARY 

DO  20  NTINF-NBOTI ,NBOTF 
20  TS(NTINF)-TINF 
DO  30  NNT-NBF,NEL 
TS(NNT)-0.0 
30  H(NNT)-0.0 

DO  40  IGT-1,NEST 
IE-ICOR( IGT,2) 

40  T(IE)-TZ 

IG-ICOR(NEST,l ) 

T( IG)-TZ 

OMEGA  IS  IN  RADIANS/HOUR 
OMEGA-RPM*2 . 0*PI*60 . 0 
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noon  non 


DO  50  KL-NB0TI,NB0TF 
50  H(KL)-HINF 
HIFN-HINF 
TSAT-TSS 
EPSO-ZOA*EZERO 
BOA*  TAL  FA 

SURFAR-NFIN* ( 2 . 0* ( S( 1 )/( 12*CALFA) ) +EPSO ) 

EPSEX-(CEXIT-(NFIN*EZERO) )/NFIN 

BETA- ( EPSEX-EPSO ) /DIV 

ZZERO- ( S ( 1 ) /12 ) /CAL FA 

AFOVAS- ( ZOA+ ( 1 . /SALFA ))/(!. +20A ) 

ZA-0 . 0 

DO  60  NSAT-1,NEST 
60  TS(NSAT)-TSAT 

TSOLID- ( TSAT+TINF )/2 . 0 

C  TEMPORARY  CHANGE  -  TFILM 

ASMOOTH-0.0 
ACASE-0.0 
QT-0 . 0 
OBJ-0.0 
QTl-0.0 
QTF-0.0 
QTRF-0.0 
QTOT-0.0 
DMTOT-0 . 0 
NR-NDIV+1 
DO  350  NZ-1,NK 

C  R  IS  THE  INCREMENTAL  CHANGE  IN  THE  RADIUS  OF  THE  CONDENSER 

R( NI ) «R2>NI *DELX*SPHI 
RB ( NI ) -RBASE>NI *DELX*SPHI 

C  EPS  IS  THE  INCREMENTAL  CHANGE  IN  THE  TROUGH  WIDTH 

EPS ( NI ) -EPSO+NI *BETA 
APS-EPS(NI) 

NODAL  POINT  COORDINATES 

CALL  COORD 
65  Z(1)-ZA 

DO  70  IZEL-1,NEFB 
NA-ICOR( IZEL,1 ) 

NB-ICOR( IZEL,2) 

XE-X(NA)-X(NB) 

YE-Y(NA)-Y(NB) 

ELZ-SQRT(XE**2+YE**2) 

70  Z( IZEL+1 )-Z( IZEL)+ELZ 

XZB-X(ICOR(NHB,l) )-X( ICOR(l,2) ) 

YZB-Y(ICOR(NHB,l) )-Y(ICOR(l,2) ) 

ZB-SQRT(XZB**2+YZB**2 ) 

ZC-ZZERO 
IM-1 

PARABOLIC  TEMPERATURE  DISTRIBUTION  ALONG  THE  FIN 
BOUNDARY, USING  LAGRANGE  INTERPOLATION 

80  TP1-T( ICOR( 1,2) ) 

TP2-T( ICOR{ NHB, 1 )  ) 

TP3-T( ICOR(NEFB, 1 ) ) 

AP1-TP1/(ZB*ZC) 

AP2-TP2/( ZB*( ZB-ZC) ) 

AP3-TP3/( ZC*( ZC-ZB) ) 
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noon  on  non  oonon 


BPl— (ZB+ZC)*AP1 
BP2— ZC*AP2 
BP3— ZB*AP3 
A1-AP1+AP2+AP3 
Bl«BPl-i-BP2-fBP3 
TC-0.0 

DO  90  NY-1, NEST 

90  TC-TC+T( IC0R(NY,2) ) 

AY- FLOAT (NY+1 ) 

TF-(TC+T( IC0R(NY,1) )+AY*TS(NY) )/(2.0*AY) 

SOLID-FLUID  PROPERTIES 
WATER  PROPERTIES 
IF(IFLUID.EQ.l)  GO  TO  91 

HPG-1093. 88-0. 5703*TS(l)+0. 00012819* (TS(1)**2) 
1-0.0000008824*(TS(1)**3) 

RHOF( NI) -62. 774-0. 00255698*TF-0.000053572*TF**2 
CF(NI )-0. 30 34+0. 0007 38927*TF-0. 000001 47 321*TF* *2 
UF(NI) -0.001397-0. 00001 4669*TF+0. 0000000631253 *TF* *2-0. 00000 
100000976569*TF**3 

FREON  PROPERTIES 

91  IF(IPLUID.EQ.O)  GO  TO  92 

HFG-69. 5459-0. 0156011*TS(l)-0.000455294*(TS(l)**2)+0. 00000104144* ( 
1TS(1)**3) 

RHOF(Nl ) -102. 059-0. 025364*TF-0. 000502649* (TF**2)+0. 00000135407 *(TF 
1**3) 

CF(NI)-0. 0594858-0. 000429765*TP+0.00000348218*TF**2-0. 000000010416 
18*TP**3 

UF(NI )-0. 00078-0. 00000525*TF+0.0000000125*TF**2 

92  UF(NI )-3600*UF(NI ) 

1F( IFIN.EQ.l)  GO  TO  93 

CW( NI) -231 .7772-0. 02222*TSOLID 

93  IF(IFIN.EQ.O)  GO  TO  94 
CW(NI) -8. 776+0. 0026 5*TSOLID 

94  CK-CW(NI) 

CONST-RHOF(NI ) **2*OMEGA**2*HFG*CPHI*CALFA*R( NI ) 

INITIAL  FILM  THICKNESS 
IF  (NI.GT.l)  GO  TO  100 

DEL(1)-1.107*( ( (TSAT-TINF)*CF(NI)/(UF(NI)*HFG) ) ** . 25 ) * ( ( UF( NI )/( 
lRHOF(NI)*OMEGA) )**0.5) 

100  CONTINUE 

AVERAGE  ELEMENT  CONVECTIVE  COEFFICIENT  ALONG 
THE  FIN  BOUNDARY 

ZSTAR-ZZERO-DEL(NI )/CALFA 
AZZ-DEL ( NI ) /SALFA 
ZZ-ZSTAR 

HDEN-( {-Al*ZZ**3)/3.0-(Bl*ZZ**2)/2.0)+ZZ*(TSAT-T( 1 ) ) 

AZS-ABS( 4*CF(NI )*UF(NI )*HDEN/CONST)**0.25 
HAC-0 .0 

DO  190  IEL-1,NEFB 
AZ-Z( lEL) 

BZ-Z( IEL+1 ) 

IF  ( ZSTAR.LE.BZ )  GO  TO  110 
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GO  TO  120 

110  IF  (HAC.NE.0.0)  GO  TO  180 
BZ-ZSTAR 

120  IF  (lEL.NE.l)  GO  TO  130 
AK-(BZ-AZ)/5.0 
ZZ-AK 
GO  TO  140 

130  AK-(BZ-AZ)/4 .0 
ZZ>*AZ 

140  ZEL-4*AK 

DO  150  NH-1,5 

HDEN-(-1.0*(Al*ZZ**3/3.0+Bl*ZZ**2/2.0) ) +ZZ* ( TSAT-T(  1 ) ) 
HZ(NH)-ABS(CF(NI)**3*C0NST/(4*0F(NI)*HDEN) )**0.25 
150  ZZ-ZZ+AK 

C  AVERAGE  H  USING  SIMPSONS  RULE 

CONH-AK*  (  HZ  ( 1 )  +  4 *HZ  (  2  )  +2 *HZ  (  3  )  +  4 *HZ  (  4  )  +HZ  (  5  )  ) /(  3 *ZEL ) 
IF  (ZSTAR.EQ.BZ)  GO  TO  160 
H( IEL)-CONH 
GO  TO  190 
160  AZ-ZSTAR 

HAZ-CONH*(AZ-Z( lEL) ) 

DELA-AZS 
170  BZ-ZdEL-t-l) 

DELB- ( BZ>ZSTAR ) *AZZ/( ZZERO-ZSTAR ) 

DELZ-  (  DELA-i>DELB )  /2 . 0 

HAC- ( BZ-AZ ) *CF ( NI ) /DELZ 

H(  IEL)-(HAZ>>>HAC)/(BZ-Z(  lEL)  ) 

GO  TO  190 
180  AZ-Z(IEL) 

DELA-DELB 
HAZ-0.0 
GO  TO  170 
190  CONTINUE 

NETI-NEFB+1 
DO  200  lEL-NETl ,NEST 
200  H( IEL)-CF(NI )/DEL(NI ) 

C 

C  ENTRY  INTO  THE  FINITE  ELEMENT  SOLUTION 

C 

CALL  FORMAF 

CALL  BANDEC  ( NSNP , NBAN , 1 ) 

C 

C  THE  TEMPERATURE  DISTRIBUTION 

C 

DO  210  NT-1, NSNP 
210  T(NT)-F(NT,1) 

TIB(NI)-T(ICOR(NBOTI,2) ) 

TT(Nl)-T(ICOR(NEFB,l) ) 

TE(NI)-T( ICOR(NEST,l) ) 

TB(NI)-T( ICOR(NBOTF, 1 ) ) 

TTS-0.0 

DO  220  NS-1,NSNP 
220  TTS-TTS+T(NS) 

PN-FLOAT(NS ) 

TSOLID-TTS/PN 

C 

C  Q  AT  THE  BOTTOM  SIDE 

C 

QBI-0.0 

DO  230  IBEL-NBOTI ,NBOTF 
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NKA«ICOR( ZBEL,1) 

NKB-IC0R(IBEL,2) 

XB-X(NKA)-X(NKB) 

YB-Y(NKA)-Y(NKB) 

ELB-SQRT(XB**2+YB**2) 

230  QBI-QBI+(T(NKA)+T(NKB)-2*TS(IBEL) ) *ELB*H( IBEL)/2.0 
QB(NI )-QBI*DELX 

ITERATION  UNTIL  CONVERGENCE  CRITERIA  IS  MET 

IF  (IM.EQ.l)  GO  TO  240 
QJ-QBI 
GO  TO  250 
240  QI-QBI 
IH-i2 

GO  TO  80 

250  AQ-ABS(QJ-QI)/QJ 

IF  (AQ.LE.CRIT)  GO  TO  260 

QI-QJ 

GO  TO  80 

260  DMDOT(NI)-2.*QBI*DELX/HFG 
DMTOT-DMTOT+DMDOT( NI ) 

Cl-RHOF(NI ) **2*OMEGA**2*R( NI ) *SPHI/( 3*UF(NI ) ) 

XCOF(l)— DMTOT 

XCOF(2)-0.0 

XCOF(3)-0.0 

XCOF(4)-Cl*EPS(NI) 

XCOF( 5)-Cl*TALFA 
H-4 

CALL  DPOLRT  (M,IER) 

IF  (ROOTR(l) .GT.O.O)  GO  TO  270 

IF  (ROOTR(2) .GT.0.0)  GO  TO  280 

IF  (ROOTR(3) .GT.0.0)  GO  TO  290 

IF  (ROOTR( 4) .GT.0.0)  GO  TO  300 

THE  CONDENSATE  THICKNESS 

270  DEL(NI+l)-ROOTR(l) 

GO  TO  310 

280  DEL(NI-»>l)«ROOTR(2) 

GO  TO  310 

290  DEL(NI+1 )-ROOTR( 3 ) 

GO  TO  310 

300  DEL(N1+1 )-ROOTR( 4 ) 

310  QEL-0.0 

IF  (NI.NE.l)  GO  TO  320 

Q  FROM  THE  TOP  SIDE 
Q  THROUGH  FIN 
QEL-0.0 

320  DO  330  IQEL-1,NEFB 
KA-ICOR( IQEL, 1 ) 

KB-ICOR( IQEL, 2 ) 

XQEL-X( KA)-X( KB) 

YQEL-Y( KB )-Y( KA) 

ELM-SQRT(XQEL**2+YQEL**2) 

QEL-QEL+( 2*TS( IQEL)-T( KA)-T(KB) )*ELM*H( IQEL)/2 . 0 
330  CONTINUE 
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QINC(NZ)-QEL*OELX 
AMTOT(NI )-DMTOT 
QET«QBL*DELX*NFIN*2 
QT-QT4-QET 
QA-QBI*DELX*NFIN*2 
C 

C  Q  THROUGH  TROUGH 

C 

QTRF-0 . 0 

DO  340  IQEL-NEFB+l.NEST 
KA-1C0R( IQEL,1) 

KB-ICOR( IQEL,2) 

XQEL-X(KA)>X(KB) 

YQBL-Y(KB)-Y(KA) 

ELM-SQRT(XQEL**2+YQEL**2) 

QTRF-QTRF+ ( 2  *TS ( 1 QEL ) -T ( KA ) -T ( KB ) ) *ELM*H ( 1 QEL ) /2 . 0 
340  CONTINUE 

QTINC ( N1 ) -QTRF*DELX 
QTOTAL(NI )-QlNC(NI )+QTINC(NI ) 

QTRFT-QTRF*DELX*NFIN*2 . 

QTF-QTF+QTRFT 

QTOT-QTOT+QA 

ASMOOTH-2  *PI *RB ( N1 ) *DELX+ASMOOTH 

ACASE-ACASE+NFIN*( ( (2*S(1) )/(12*CALFA) )+EPSO)*DELX 
350  CONTINUE 

C  EVALUATE  OBJECTIVE  FUNCTIONS  AND  CONSTRAINTS 

OBJ— (QTOT) 

WRITE (15,*)  'OBJECTIVE  FUNCTION-',  OBJ,'BFIN-',  S(l) 

C  FIRST  CONSTRAINT  IS  TO  ENSURE  THE  RATIO  ZOA  IS  NOT  NEGATIVE 

G(l)— ( ( (CBASE-(2*NFIN*S(1)*TALFA/12) )/NFIN)/(S{l)*TALFA/6) ) 
G(1)-1000.0*G(1) 

C  THE  SECOND  CONSTRAINT  IS  TO  ENSURE  THE  CONDENSATE  LEVEL  IS  NOT 

C  GREATER  THAN  THE  FIN  HEIGHT 
G(2)— ( (S(1)/12.0)-DEL(NI) ) 

XPLOT(NO)-S(l) 

FNOB  J  ( NO )  —OBJ 

ARAT 1 0- ACAS  E/AS MOOTH 

WRITE(13,*)  XPLOT(NO) ,FNOBJ(NO) 

WRITE(14,*)  ARATIO,  FNOBJ(NO) 

GO  TO  10 
360  CONTINUE 
C 

BFIN-S(l) 

C  *****  OUTPUT  NODE  ***** 

WRITE  (15,630) 

DO  370  NR-1,NK 

370  WRITE  (15,640)  NR,QINC( NR ) ,QTINC( NR ) , QTOTAL ( NR ) 

WRITE  (15,650)  QT,QTF 

WRITE  ( 15,660)  CLI , CANGL ,RBASEI , R2l , THICKI , BFIN , RPM, TSS , TINF , HINF , 
IRIT, FANGL,ZOA, IFF 
WRITE  (15,661)  AFOVAS 

WRITE  (15,670)  BOA,ZOA,NFIN,BVIN,SURFAR 
WRITE  (15,680) 

DO  380  NP-1,NSNP 
TCC(NP)-. 5555555* (T(NP) -32 ) 

380  WRITE  (15,690)  NP , X( NP ) , Y( NP ) , T( NP ) , TCC( NP ) 

WRITE  (15,700) 

DO  390  KKL-l,NBOTF 
NKX-ICOR(KKL,l) 

NKY-ICOR(KKL,2) 
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YP-Y(NKX)-Y(NKY) 

EXY-SQRT( XP**2+YP**2 ) 

QEP-ABS ( ( T( NKX ) +T( NKY ) -2*TS ( KKL ) ) *EXY*H ( KKL)/2 . 0 ) 

QEP-QEP*DELX 

390  WRITE  (15,710)  KKL , H( KKL) , EXY,QEP 
WRITE  (15,720)  CRIT 

WRITE  (15,730)  HFG,NFIN,H(NBOTF) ,TSAT,RPM,QTOT.QT,FANGL 
WRITE( 15,734 ) 

WRITE(15,735)  ROOTR(l) ,ROOTI(l) 

WRITE(15,735)  ROOTR( 2 ) , ROOTI ( 2 ) 

WRITE(15,735)  ROOTR( 3 ) , ROOTI ( 3 ) 

WRITE(15,735)  ROOTR( 4 ) , ROOTI ( 4 ) 

WRITE  (15,740) 

DO  400  NR-1,NK 

400  WRITE  (15,750)  NR,DEL(NR) ,QB(NR) ,AMTOT(NR) ,TIB(NR) ,TT(NR) ,TE(NR) , 
1TB (NR) 

WRITE  (15,760) 

LO  410  NG-1,NDIV,2 

410  WRITE  (15,770)  NG,G4(NG) ,CF(NG) ,RHOF(NG) ,UF(NG) ,EPS(NG) ,R(NG) , 
IQINC(NG) 

RETURN 

412  FORMAT  ( 8X, E12 . 5 , 8X, E12 . 5 ) 

420  FORMAT  (5l5) 

430  FORMAT  (/2X, 15HNO. OF . ELEMENTS-, I5 , lOX, 18HNO .OF . SYSTEM  N.P.-, 
1I5,/,2X,13HN0.0F  BANDED-, 15) 

435  FORMAT  (/2X, ' IFLUID-' , 15 , lOX, ' IFIN-S I5 ) 

436  FORMAT  (2X,'IFLUID  -  0  FOR  WATER,  AND  1  FOR  FREON') 

437  FORMAT  (2X,'IFIN  -  0  FOR  COPPER,  AND  1  FOR  STAINLESS  STEEL') 

440  FORMAT  (4I5) 

450  FORMAT  ( /2X, 7HELEMENT, lOX, 3HNP1 , 14X, 3HNP2 , 15X, 3HNP3 ) 

460  FORMAT  (7G10.5) 

470  FORMAT  (4X,5HCLI-  , E12 . 5 ,/, 4X, 7HCANGL-  , E12 . 5 ,/, 4X, 8HRBASEI . - , E12 . 
15,/,4X,5HR2I-  ,E12.5,/,4X,8HTHICKI-  , E12 . 5 ,/, 4X, 6HBFIN-  ,El2.5,/,4 
2X,4HTZ-  ,E12.5) 

480  FORMAT  (5l5) 

490  FORMAT  (4X,6HNDIV-  , IlO ,/, 4X, 6HNEST-  , IlO ,/, 4X, 6HNEFB-  ,I10,/,4X,7 
IHNBOTI-  ,I10,/,4X,7HNBOTF-  ,110) 

500  FORMAT  (4F10.2) 

510  FORMAT  (4X,5HRPM-  , E12 . 5 ,/, 4X, 5HTSS-  , El2 . 5 ,/, 4X, 6HTINF-  ,E12.5,/, 
14X,6HHINF-  ,E12.5) 

520  FORMAT  (G10.9) 

530  FORMAT  (4X,6HCRIT-  ,E12.5) 

540  FORMAT  (2G10.5) 

550  FORMAT  (4X,7HFANGL-  , E12 . 5 ,/, 4X, 6HNFIN-  ,15) 

560  FORMAT  (I5) 

570  FORMAT  (4X,5HIFF-  ,110) 

580  FORMAT  (16l5) 

590  FORMAT  (6l5) 

600  FORMAT  ( ///5X , 4HTIB- , 1 5 , lOX, 3HTT- , 1 5 ,/, 5X , 3HTE- , 1 5 , / 

1,6X,3HTB-,I5) 

610  FORMAT  ( //I OX , 17HCRASH , CRASH , CRASH ) 

620  FORMAT  (//5X, 4 ( E12 . 7 , 3X) ) 

630  FORMAT  ( 2X , 7HSTATION , 2X, 4HQFIN , 17X , 7HQTROUGH , 15X , 6HQTOTAL ) 

640  FORMAT  ( 4X , 1 5 , E12 . 5 , 1 OX , E12 . 5 , lOX , E12 . 5 ) 

650  FORMAT  (//, 4X, llHQFIN  TOTAL-, El2 . 5 , lOX, 15HQTROUGH  TOTAL-  ,E12.5) 
660  FORMAT  (/////, 4X , 5HCLI-  , E12 . 5 , 5X, 7HCANGL-  , E12 . 5 ,/, 4X, 8HRBASEI-  , 
1E12.5,2X,5HR2I-  , E12 . 5 ,/, 4X, 8HTHICKI-  , E12 . 5 , 2X , 6HBFIN-  ,El2.5,/,4 
2X,5HRPM-  ,E12.5,5X,5HTSS-  , E12 . 5 ,/, 4X, 6HTINF-  , E12 . 5 , 4X , 6HHINF-  ,E 


312.5,/,4X,6HCRIT-  , E12 . 5 , 4X, 7HFANGL-  , E12 . 5 ,/, 4X, 5HZ0A-  ,E12.5,5X, 
45HIFP-  ,110) 

661  FORMAT(4X, 'FIN  AREA/SMOOTH  AREA-' , El2 . 5 ) 

670  FORMAT  ( IHI ,//2X , 4HBOA- , G12 . 5 , 5X, 4HZOA- , G12 . 5 , 3X , 5HNFIN- , 1 5 , 5X , 

1/, 5HBVIN- , G12 . 5 , 5X, 1 3HSURFACE  AREA- , G12 . 5 ) 

680  FORMAT  ( //5X, 2HNP , 6X, IHX, 12X, IHY , 12X, IHT , 12X, 2HTC ) 

690  FORMAT  ( /2X, 1 3 , 3X , 4 ( FlO . 6 , 3X ) ) 

700  FORMAT  ( //2X , 2HEL , 8X , IHH , llX , 9HEL-LENGTH , 1 5X , 4HQ-EL) 

710  FORMAT  (/2X, 12 , 3X, E12 . 5 , 3X, E12 . 5 , lOX, E12 . 5 ) 

720  FORMAT  ( /2X, 22HCONVERGENCE  CRITERIAN- , El 5 . 8 ) 

730  FORMAT  ( IH  , // , 5X , 4HHFG- , El 2 . 5 ,/, 5X , 1 IHNO . OF  FINS- , I  5 , / , 5X , 

1 6HH-OUT- , E12 . 5 , / , 5X , 5HTSAT- , E12 . 5 , / , 5X , 4HRPM- , El 2 . 5 , / , 5X , 6HQ-BOT= , 
2E12 . 5 ,/, 5X, 6HQFIN  E12 . 5 ,/, 5X, llHHALF-ANGLE- , F8 . 3 ) 

734  FORMAT (/5X, 'ROOTS: ' ,5X, 'REAL  PARTS ', 1 5X ,' IMAGINARY  PARTS') 

735  F0RMAT(15X,E12.5,15X,E12.5) 

740  FORMAT  ( IHO , 6X, IHJ , 4X, 14HFILM  THICKNESS , 8X, 2HQB , lOX, 8HMASS-TOT , 
1/,4X, 3HTIB,8X,2HTT,10X,2HTE,8X,2HTB) 

750  FORMAT  ( IH  , 4X, 14 , 4X, F12 . 10 , 4X, FlO . 4 , 6X , F9 . 5 , 6X, F5 . 1 , /, 
16X,F5.1,6X,F5.1,6X,F5.1) 

760  FORMAT  ( IHO , 6X,1HJ,6X,6HK-WALL,4X,6HK-FILM,3X,7HDENSITY,4X,9HVISC- 
IFILM,/, 6X, 7HEPSILON, 5X, 6HRADIUS, 5X, 4HQINC) 

770  FORMAT  ( IH  , 4X, 14 , 4X, F7 . 3 , 4X, F6 . 4 , 4X, F6 . 3 , 4X, F9 . 7 , 

1/, 4X , F9 . 7 , 4X , F7 . 5 , 5X , F5 . 1 , IX , F7 . 3 ) 

END 


THIS  SUBROUTINE  DEFINES  THE  POSITIONS  OF  SYSTEM  COORDINATE  POINTS 

SUBROUTINE  COORD 

COMMON/ADS/DF( 21 ) ,G( 10 ) , IDG ( 100 ) , IGRAD, INFO, lOPT, lONED, I PRINT, 
:ISTRAT,IWK(2000) ,IZ( 30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W(21, 30) ,WK( 5000) , 
:  NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A( 200,50) , AMTOT( 200 ) , APS , B ( 3 ) , BFIN, BOA, BVIN, C( 3 ) , 
:CANGL,CF( 200) , CK , CLI , COF ( 5 ) , CW( 200) , DEL ( 200 ) ,DMDOT( 200),EA(3,3), 

:  EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ( 100 ) ,H(200) , HINF, HZ (  200  ) , 

:QB( 200) ,QINC( 200) ,QTINC( 200 ) ,QTOT,QTOTAL( 100 ) , R( 200 ) , RB( 200 ) , 
;RBASEI ,R2I,RHOF( 200) ,ROOTI( 4) ,RPM,ROOTR( 4 ) ,T( 200 ) , TALFA, TB ( 200 ) , 
:TCC(200) ,TE(200) , THICK, THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS, 
:TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5),XPLOT(200) , Y( 200 ) , Z ( 200 ) ,ZOA, 
:DOBF,DOTH, ICOR( 200, 3 ) , IFF, JINT, JLC, JTC, KFF( 50 ) ,KFIN( 50 ) ,KT,NBAN, 
:NEL,NFIN,NSNP 


C  DELH  IS  THE  STANDARD  DIVISION  OF  FIN  HEIGHT 

DELH-S( 1 )/( 12*DOBF) 

X( 1 )-0.0 

Y( 1 )-THICK+( S(1 )/12 ) 

N-1 

DO  20  I-1,IFF 
ICA-KFIN( I ) 

ICB-KFF( I ) 

CBA- FLOAT ( ICB-ICA) 

AN-0 . 0 

DO  10  II-ICA,ICB 
X( II )-X( 1 )+N*AN*DELH*TALFA/CBA 
Y( II )-Y( 1 )-N*DELH 
10  AN-AN+1.0 
20  N-N+1 
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AN-0.0 

ICD-ICB-ICA+1 

DO  50  J-JTC, JLC, JINT 

X(J)-X(1) 

Y( J)-( 1 .0-AN/DOTH)*THICK 
DO  30  JJ-1,ICD 

X( J+JJ )-X( J ) +JJ*EZERO/( 2* ( CBA+1 . 0 ) ) 
30  Y(J+JJ)-Y(J) 

JJ-ICD 

DO  40  K-1,KT 

X( J+JJ+K)-X( J+JJ ) +K*APS/( 2 . 0*KT) 

40  Y( J+JJ+K)-Y( J) 

50  AN-AN+1.0 
RETURN 
END 


*  THIS  SUBROUTINUE  IS  USED  TO  FORMULATE  THE  FEM  EQUATIONS 

* 

SUBROUTINE  FORMAF 

* 

COMMON/ADS/DF(21) ,G(10) ,IDG(100) ,IGRAD,INFO,IOPT,IONED,IPRINT, 
;ISTRAT,IWK( 2000) ,12(30) ,OBJ,S( 2 ) ,VLB( 2 ) ,VUB( 2) , W( 21 , 30 ) ,WK( 5000 ) 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

* 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B(3) ,BFIN,BOA,BVIN,C{3) , 
:CANGL,CF(200) ,CK,CLI,COF(5) ,CW(200) ,DEL(200) ,ONDOT(200) ,EA(3,3) , 
: EPS ( 200) ,EZERO,F( 200,1) ,FANGL,FNOBJ( 100) ,H(200) ,HINF, HZ (200) , 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL(100) ,R(200) ,RB(200) , 
;RBASEI,R2I,RHOF(200) ,ROOTI(4) ,RPM,ROOTR(4) ,T(200) ,TALFA,TB(200) , 
;TCC(200) ,TE( 200) , THICK, THICKI, TIB (200) ,TINP,TS( 200) ,TSAT,TSS,. 
:TT(200) ,T2,UF(200) ,X(200) ,XCOr(5),XPLOT(200) ,Y(200) ,Z(200) ,20A, 
;DOBF,DOTH,ICOR(200,3) ,IFF,JINT,JLC,JTC,KPF(50) ,KFIN(50) ,KT,NBAN, 
:NEL,NFIN,NSNP 


DO  20  N-il,NSNP 
F(N,l)-0.0 
DO  10  HA-1, NBAN 
10  A(N,nA)-0.0 
20  CONTINUE 

DO  70  IEL-1,NEL 
.  IA-ICOR( IEL,1) 

IB-ICOR( IEL,2) 

IC-ICOR( lEL, 3  ) 

B(1)-Y(IB)-Y(IC) 

B(2)-Y( IC)-Y( lA) 

B( 3)-Y( IA)-Y( IB) 

C(1)-X(IC)-X(IB) 

C(2)-X( IA)-X( IC) 

C( 3 )-X( IB)-X( lA) 

C  LENGTH  BETWEEN  ELEMENT  NODES  1  AND  2 

EL-SQRT(C( 3)**2+B( 3)**2) 

C  AREA  OF  TRIANGULAR  ELEMENT 

AS-ABS( (B(1)*C(2)-B(2)*C(1) )/2.0) 
HC-H( IEL)/CK 
DO  60  J-1,3 
JJ-ICOR( lEL, J) 

DO  50  K-1,3 
KK-ICOR( IEL,K) 
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C  FORMING  THE  A  MATRIX 

EA( J,K)-(B( J)*B(K)+C( J)*C(K) )/(4*AS) 

IF  (HC.EQ.0.0)  GO  TO  40 

HEL-HC*EL/6.0 

IF  (J.EQ.3)  GO  TO  40 

IF  (K.EQ.3)  GO  TO  40 

IF  (J.EQ.K)  GO  TO  30 

EA( J,K)-EA( J,K)+HEL 

GO  TO  40 

30  EA( J,K)-EA( J,K)+2*HEL 
40  IF  (KK.LT.JJ)  GO  TO  50 
NW-KK-JJ+1 

A( JJ,NW)-A( JJ,NW)+EA( J,K) 

50  CONTINUE 
60  CONTINUE 

C  FORMING  THE  F  MATRIX 

FE-HC*TS ( lEL) *EL/2 . 0 
F(IA,1)-F(IA,1)+FE 
F(IB,1)-F(IB,1)+FE 
70  CONTINUE 
RETURN 
END 


*  EQUATION  SOLVER  OF  A  SYMMETRIC  MATRIX  THAT  HAS  BEEN  TRANS- 

*  FORMED  INTO  BANDED  FORM. 

* 

SUBROUTINE  BANDEC  ( NEQ , MAXB , NVEC ) 

* 

COMMON/ADS/DF( 21 ) ,G( 10 ) , IDG( 100) , IGRAD, INFO, lOPT, lONED, IPRINT, 
:ISTRAT,IWK(2000) ,IZ(30) ,OBJ,S(2),VLB(2) ,VUB(2) ,W(21,30) ,WK(5000) , 
; NCOLA , NCON , NDV , NGT , NRA , NRl WK , NRWK 

* 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B(3) ,BFIN,BOA,BVIN,C( 3) , 
:CANGL,CF(200) , CK , CLI , COF( 5 ) , CW( 200 ) , DEL ( 200 ) , DMDOT( 200 ) , EA( 3 , 3 ) , 

: EPS (200) ,EZERO,F( 200,1) , FANGL , FNOB J ( 1 0 0 ) ,H(200) , HINF , HZ  (  200 ) , 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL(100) , R( 200 ) , RB ( 200 ) , 
:RBASEI,R2l,RHOF( 200) , ROOTI ( 4 ) , RPM, ROOTR( 4 ) , T( 200) , TALFA, TB ( 200 ) , 
:TCC(200) ,TE( 200) , THICK, THICKI, TIB (200) ,TINF,TS( 200) ,TSAT,TSS, 
:TT(200) ,TZ,UF(200) ,X( 200 ) ,XCOF( 5 ) ,XPLOT( 200 ) , Y( 200 ) ,Z(200) ,ZOA, 
:DOBP,DOTH,ICOR( 200,3) , IFF , JINT, JLC, JTC, KFF( 50) ,KFIN( 50 ) ,KT,NBAN, 
:NEL,NFIN,NSNP 


LOOP-NEQ-1 
DO  20  I-l,LOOP 
MB-I+1 

NB-MINO ( I+MAXB-1 , NEQ ) 

DO  20  J-NB,NB 

L-J+2-MB 

D-A(I,L)/A(I,1) 

DO  10  MM-1,NVEC 
10  F(  J,MM)-F(  J,MM)-D’»F(  I  ,MM) 
MM-MINO ( MAXB-L+1 , NEQ- J+1 ) 

DO  20  K-1,MM 
NN-L+K-1 

20  A( J,K)«A( J,K)-D*A( I,NN) 

DO  30  I-1,NVEC 

30  F(NEQ, I )-F(NEQ, I )/A(NEQ,l ) 
DO  50  1-2, NEQ 
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J-NEQ-I+1 

K-MINO ( NEQ-J+1 , MAXB ) 

DO  50  MM-1,NVEC 
DO  40  L-2,K 
MB-J+L-1 

40  F( J,MM)-F( J,MM)-A( J,L) *F(MB,MM) 
50  F( J,MM)-F( J,MM)/A( J,l) 

RETURN 

END 


SUBROUTINE  DPOLRT  COMPUTES  THE  ROOTS  OF  A  REAL 
POLYNOMIAL  USING  A  NEWTON-RAPHSON  ITERATIVE 
TECHNIQUE. 


SUBROUTINE  DPOLRT  (H,IER) 

COMMON/ADS/DF( 21 ) ,G(10) ,IDG(100) , IGRAD, INFO, lOPT, lONED, IPRINT, 
:ISTRAT,IWK(2000),I2(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W( 21 , 30 ) ,WK( 5000 ) , 
;  NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A( 200,50) ,AMTOT(200) ,APS,B( 3) , BFIN , BOA , BVIN , C ( 3 ) , 
:CANGL,CF(200) , CK, CLl , COF ( 5 ) , CW( 200 ) , DEL (200) ,DMDOT( 200 ) ,EA( 3, 3) , 
:EPS ( 200 ) ,EZERO,F( 200,1 ) ,FANGL,FNOBJ( 100) ,H( 200) , HINF , HZ ( 200 ) , 
;QB(200) ,QINC(200) ,QTINC(200),QTOT,QTOTAL(100) ,R(200) ,RB(200) , 
:RBASEI,R2I,RHOF(200) ,ROOTl(4) ,RPH,ROOTR( 4) ,T( 200 ) ,TALFA,TB( 200 ) , 
:TCC(200)  ,TE(  200)  , THICK, TH1CKI,T1B(  200.)  ,T1NF,TS(  200)  ,TSAT,TSS, 
:TT(200) ,T2,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200) ,Z(200) ,ZOA, 
;DOBF,DOTH,ICOR(200,3) ,IFF,JINT, JLC, JTC,KFF(50) ,KF1N( 50) ,KT,NBAN, 
:NEL,NFZN,NSNP 


lFlT-0 

N-M 

IER-0 

IF  (XCOF(N+l))  10,40,10 
10  IF  (N)  20,20,60 

SET  ERROR  CODE  TO  1 

20  IER-1 
30  RETURN 

SET  ERROR  CODE  TO  4 

40  IER«4 
GO  TO  30 

SET  ERROR  CODE  TO  2 

50  IER-2 
GO  TO  30 

60  IF  (N-36)  70,70,50 
70  NX-N 
NXX-N+1 
N2-1 
KJl-N+1 
DO  80  L-1,KJ1 
MT-KJl-L+1 
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80  C0F(MT)-XC0F(L) 

SET  INITIAL  VALUES 

90  X0-. 00500101 
YO-0. 01000101 

ZERO  INITIAL  VALUE  COUNTER 
IN-0 

100  XX-XO 

INCREMENT  INITIAL  VALUES  AND  COUNTER 

XO— 10.0*YO 
YO— 10 .0*XX 

SET  X  AND  Y  TO  CURRENT  VALUE 

XX-XO 
YY-YO 
IN-IN+1 
GO  TO  120 
110  IFIT-1 
XPR-XX 
YPR-YY 

EVALUATE  POLYNOMIAL  AND  DERIVATIVES 

120  ICT-0 
130  UX-0.0  . 

UY-0 . 0 
/-O.O 
YT-0.0 
XT-1.0 
U-COF(N+l ) 

IF  (U)  140,270,140 
140  DO  150  I-1,N 
L-N-I+1 

XT2-XX*XT-YY*YT 
YT2-XX*YT+YY*XT 
U-U+COF(L)*XT2 
V-V+COF( L) *YT2 
FI-I 

UX-UX+FI*XT*COF( L) 

UY-UY-FI * YT*COF ( L ) 

XT-XT2 
150  YT-YT2 

SUMSQ-UX*UX+UY*UY 
IF  (SUMSQ)  160,230,160 
160  DX-(V*UY-U*UX) /SUMSQ 
XX-XX+DX 

DY—  (  U*UY+V*UX  )  /SUMSQ 
YY-YY+DY 

IF  (ABS(DY)+ABS(DX)-1 .OE-05)  210,170,170 

STEP  ITERATION  COUNTER 

170  ICT-ICT+1 

IF  (ICT-500)  130,180,180 
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180  IF  (IFIT)  210,190,210 
190  IF  (IN-5)  100,200,200 

SET  ERROR  CODE  TO  3 

200  IER-3 

GO  TO  30 

210  DO  220  L-1,NXX 
MT-KJl-L+1 
TEMP-XCOF(MT) 

XCOF(MT)-COF(L) 

220  COF(L)-TEMP 
ITEMP-N 
N-NX 
NX- I TEMP 

IF  (IFIT)  250,110,250 
230  IF  (IFIT)  240,100,240 
240  XX-XPR 
YY-YPR 
250  IFIT-0 

IF  (ABS(YY/XX)-1.0E-04)  280,260,260 
260  ALPHA- XX-i- XX 

SUMSQ-XX*XX+YY*YY 
N-N-2 
GO  TO  290 
270  XX-0.0 
NX-NX-1 
NXX-NXX-1 
280  YY-0.0 

SUHSQ-0.0 

ALPHA-XX 

N-N-1 

290  COP(2)-COF(2)+ALPHA*COF(l) 

DO  300  L-2,N 

300  COF( L+1 )-COF( L+1 )+ALPHA*COF( L)-SUMSQ*COF( L-1 ) 
310  ROOTI(N2)-YY 
ROOTR(N2)-XX 
N2-N2-fl 

IF  (SUHSQ)  320,330,320 
320  YY— YY 

SUMSQ-0.0 
GO  TO  310 

330  IF  (N)  30,30,90 
END 
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